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We examine whether the cosmic ray positron excess observed by PAMELA can be explained 
by neutrahno annihilation in the Next-to-Minimal Supersymmetric Standard Model (NMSSM). 
The main dark matter annihilation products are the lightest CP-even scalar hi plus the lightest 
CP-odd scalar ai, with the ai decaying into two muons. The energetic positrons needed to explain 
PAMELA are thus obtained in the NMSSM simply from kinematics. The required large annihilation 
cross section is obtained from an s-channel resonance with the heavier CP-odd scalar 02. Various 
experiments constrain the PAMELA-favored NMSSM parameter space, including collider searches 
for a light ai. These constraints point to a unique corner of the NMSSM parameter space, having 
a lightest neutralino mass around 160 GeV and a very light pseudoscalar mass less than a GeV. 
A simple parameterized formula for the charge-dependent solar modulation effects reconciles the 
discrepancy between the PAMELA data and the estimated background at lower energies. We 
also discuss the electron and gamma ray spectra from the Fermi LAT observations, and point out 
the discrepancy between the NMSSM predictions and Fermi LAT preliminary results and possible 
resolution. An NMSSM explanation of PAMELA makes three striking and uniquely correlated 
predictions: the rise in the PAMELA positron spectrum will turn over at around 70 GeV, the dark 
matter particle mass is less than the top quark mass, and a light sub-GeV pseudoscalar will be 
discovered at colliders. 
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I. INTRODUCTION 

Recently, the PAMELA collaboration has ob- 
served an anomalous positron abundance in cos- 
mic radiation The positron over electron frac- 
tion turns over and appears to rise at energies from 
10 GeV to 100 GeV. However, from the same de- 
tector, no obvious antiproton excess is seen for the 
same energy range Many suggestions have been 
made to explain the positron excess at PAMELA. 
Among different approaches, dark matter annihila- 
tion is especially interesting and could imply future 
signals in dark matter direct detection experiments 
and/or at the Large Hadron Collider (LHC). 

In recent model-independent studies, dark matter 
particles are required both to annihilate dominantly 
to leptons and to have a much larger annihilation 
rate in the galactic halo than would be implied by 
a traditional thermal relic estimate Q j^. The lack 
of antiproton excess in the PAMELA experiment [1] 



can be explained if dark matter particles first an- 
nihilate into some intermediate particles, which are 
so light that their decays to hadrons are kinemat- 
ically forbidden 0], or couple dominantly to Stan- 
dard Model (SM) leptons Q. To explain the large 
dark matter annihilation cross section in the galac- 
tic halo while remaining consistent with a thermal 
dark matter relic abundance, one can introduce an 
attractive force between two dark matter particles 
and use the Sommerfeld enhancement to boost the 
annihilation cross section in the galactic halo 
This scenario suggests interesting signatures [1] and 
can be tested at the LHC ^9^ . A second approach is 
to consider non-thermal relics. Other long-lived par- 
ticles can decay to the dark matter particles and in- 
crease the dark matter relic abundance in the galac- 
tic halo 10]. 

Instead of constructing a dark matter model ad 
hoc to explain the cosmic ray observations, in this 
paper we ask whether there is an existing well- 
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motivated model that naturally contains the nec- 
essary ingredients to explain the positron excess of 
PAMELA. The most developed framework to ad- 
dress the naturalness problem of the SM is the 
Minimal Supersymmetric Standard Model (MSSM), 
which provides the lightest superpartner (LSP) pro- 
tected by R-parity as the dark matter candidate. 
In general, the annihilation products in the MSSM 
contain not only leptons but also a large fraction of 
hadrons. This makes it difficult for the MSSM to 
explain the positron excess at PAMELA, and at the 
same time be consistent with the antiproton spec- 
trum at PAMELA [111- From the theoretic side, the 
MSSM suffers the /i-problem, which can be solved 
elegantly by introducing a new gauge singlet chiral 
supermultiplet, as proposed in the Next-to-Minimal 
Supersymmetric Standard Model [12]. A recent ex- 
ploration of the NMSSM parameter space shows that 
the lightest CP-odd particle ai (mainly from the sin- 
glet component) can be naturally^ lighter than 2 mi, 
and mainly decay into two r's [13]. One should no- 
tice that the mass of ai in the NMSSM is protected 
by the f/(l)fl symmetry and can even be lighter than 
1 GeV if the soft terms associated with the singlet 
are small. In this case the ai will decay mainly into 
two muons for a mass of a few hundred MeV. 

Therefore, if the dark matter candidate in the 
NMSSM can annihilate mostly into ai's, we can have 
a leptonic final state in the annihilation products 
simply from kinematics. The NMSSM provides the 
necessary ingredients to make this happen. Notice 
that the dark matter candidate in the NMSSM is 
the LSP neutralino (for existing studies for light 
neutrahno dark matter in the NMSSM, see [l^). 
If its mass happens to be around half that of the 
heavier CP-odd scalar 02 mass, a large dark mat- 
ter annihilation cross section is obtained through 
the s-channel resonance effects with 02. To have 
leptons dominant in the final state, we should have 
a large branching ratio of 02 to ai plus hi. This 
will naturally happen, provided that the dark mat- 
ter LSP mass is less than the top quark mass, and 
thus that the decay of 02 to tt is kinematically for- 
bidden. Since the CP-odd scalar coupling to other 
fermions is proportional to their Yukawa couplings. 



the final state can dominantly be ai -I- hi with the 
former decaying to leptons. The kinematics helps us 
to obtain hard leptons in the final state of the dark 
matter annihilations. 

In Section [Til we develop the notation by deriv- 
ing the spectrum and interactions in the NMSSM, 
and show two sets of representative model points 
allowed by current experimental constraints. We 
calculate the dark matter annihilation cross section 
in Section lllll and positron excesses from neutralino 
annihilation in Section IIVI In Section |Vl we con- 
sider the constraints on the model parameter space 
from the PAMELA antiproton spectrum. We discuss 
the gamma-ray spectrum in Section IVII and point 
out a discrepancy with recent preliminary Fermi 
LAT results and comment on a possible resolution. 
In Section IVIII we demonstrate the consistency of 
our model points with direct constraints from LEP, 
Tevatron, CLEO, B-factories and the magnetic mo- 
ment of the muon. Finally, we discuss dark matter 
direct detection and conclude in Section IVlIII 



II. SPECTRUM AND INTERACTIONS IN 
THE NMSSM 

To describe the NMSSM model, we follow the no- 
tation in the Ref. [l^. The superpotential in the 
NMSSM is 

W ^ XSHuHd + '^S\ (1) 

and the soft supersymmetry-breaking terms are 

V = XA^SH^Hd + ^A,.S^ + h.c. . (2) 

Here, hatted capital letters denote superfields, and 
unhatted capital letters the corresponding scalar 
components. The minimization of the scalar po- 
tential determines their vacuum expectation values 
(VEVs): h^ = (H^), hd = (Hd) and s = {S). The 
electroweak scale v = ^/h'^ + h"^ — 174 GeV. Since 
fi^g = A s, there are four new parameters in the 
NMSSM, which we take to be real: A, A\, k and 
A^. With sign conventions for the fields, A and 
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tan/3 = hu/hd are positive, while Ax, k, and 
^eff can have either sign. 

There exists a symmetry for the NMSSM, 
which is spontaneously broken and induces a domain 
wall problem. One can introduce higher dimension 
operators to explicitly break this discrete symme- 
try and perhaps circumvent this problem . Since 
those Z3 breaking operators have small effects on 
the analysis performed in this paper, we will neglect 
them from now on. 



A. Neutralinos 

There are five neutralinos in the NMSSM: the 
U{\)y gaugino Ai, the neutral SU{2)w gaugino A2, 
the Higgsinos V'S and V'S and the singlino -0s. In 
the basis — (— «Ai, — jA2, V'u: V'dj "^s), we have the 
neutralino mass matrix 



where 
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32 



-\hd 



fll hd 
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92 hd 
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-Xhd 

2 KS 



•(4) 



To obtain the needed dark matter annihilation rate, 
we consider the case that the LSP, the lightest neu- 
tralino, is mainly made of the bino with mixings 
with Higgsinos. For gtV,Xv < |Afi| < |^| < IM2I, 
a moderate tan (3 > 1 and /i < 0, the lightest neu- 
tralino is approximately 

X ~ — iAi — (sinai cos/3 — sin 012 sin/3)-0° 

— (sinQ;2 cos/3 + sinai sin/3)';/'d, (5) 



with 
ai 



— arctan — 



V2gi V fi cos 2/3 



— arctan 
2 



V2giv{Mi 



(A2 - gf/2)v^ 
+ ji sin 2/3) 



(6) 
(7) 



and with its mass approximated as 

gf w^(Mi + fi sin 2/3) 



TTiy = All 



2 {Ml - ^^^) 



(8) 



Hereafter, we use a simple notation x to replace the 
usual notation xl for the lightest neutralino. 



B. Higgs sector at tree level 

The charged Higgs = cos f3 + sin /3 
has a mass 



sin 2/3 

For n and A of order of unity, the mass oi is 
generically 0[^jl). 

Expanding around the Higgs fields VEVs, the neu- 
tral scalar fields are defined as 



Hu — 



Hd — hd 



S 



uR 



V2 

HdR + i Hdi 



V2 

Sb + i Si 
V2 



(10) 



For the three CP-even neutral states, we can diag- 
onalize their 3x3 mass matrix by an orthogonal 
matrix Sij to obtain the mass eigenstates (ordered 
in mass): hi = Sij {Hur, HdR, Sii)j , with masses 
denoted by to/i^ . For w <C s, we write the lightest 
CP-even Higgs as 

hi « cos a cos 9s HuR + sin a cos 9s HdR. 

-sm0sSR, (11) 

with 

TT ^ A sin 2/3 sin 4/3 

- - ^ (1 - ^ sin2/3) + 0{vys'){12) 

The singlet component of hi is small and suppressed 
by v/s. 
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There are three CP-odd pseudoscalar fields, one 
of which is a massless Goldstone mode eaten by the 
Z boson. Dropping the Goldstone mode, the re- 



maining 2x2 mass matrix in the (A, Sj) basis with 
A = cos /3 Huj + sin /3 Hdi , is 



^odd 



+ K{A^ -2ks) 4XKhuhd + XAx^ - SkA^s 



(13) 



We introduce a mixing angle 6a to diagonalize the above matrix: 

2s{Ax + Ks) 2 /, v'^ {Ax - 2k sin^ 2/3 



tan 9 a 



V {Ax - 2ks) sin 2/3 



cos 0A 



As'^{Ax + Ksy + v'^{Ax - 2ks)2 sm'^2l3 ' 
I 



(14) 



and arrive at the physical CP-odd states (ordered 
in mass) 

Oi = cos 9 a{cos l3 H^i + sin/3ffd/) + sin 0^5/, 
a2 = — sin^A(cos/3-ff„/ -|- sin/3-H'd/) + cosOaSj. 

(15) 

The lightest CP-odd particle ai is mainly composed 
of the singlet field when cos 6a ^ 0, and is a doublet 
field otherwise. Since we are looking for a very light 
scalar, wc observe from the determinant of Alodd 
that this occurs if and Ax are small. This is 
technically natural since ^ is a symmetry- 

enhancing limit. 

The heavier CP-odd scalar mass is approximately 



and 



2As(^a + Ks) 



+ 



sin 2/3 

Xv'^ {Ax - 2ksY sin 2/3 



(16) 



2 s {Ax + ks) 

with the condition Ax + ks > 0. 

For small Ax and A^, we have the following ap- 
proximate formulae: 

3s(3AAAf2 sin2/3 - 2kA,,s'^) 



Mi 



Mi 



2kXs^ 
sin 2/3 



2(s2 -I- t;2 sin^ 2/3) 
-I- 2kXv'^ sin2/3. 



(17) 
(18) 



tan^^ 



cos^ 6 a 



V sin 2/3 

sin^ 2 P 
sin^ 2 /3 -t- s2 



w2 sin^ 2 /3 



(19) 



for V <^ s. 



C. Interactions and decay modes 

Since we are interested in dark matter annihilation 
through s-channel 02 exchange, we list the relevant 

vertices associated with CP-odd scalars in this sec- 
tion. For V <^ s, the couplings of CP-odd scalars to 
fermions are 



(I2 tL tji 
a2 bL b% 

dl tL 

ai bL b% 



— I 



. rrit sin^^ 



V2v tan/3 ' 
mb tan /3 sin 6a 

. mt cos 6 A 



\f2v tan/3 ' 

m(, tan/3 cos 6^ 
* F 



(20) 
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The couplings of 02 and ai to the hghtest neutraUno 
are 

i ga2xx = *5i sin 0^ (cos 2/3 sinai — sin 2/3 sina2), 

(21) 

i gaixx = cos (cos 2/3 sinai — sin 2/3 sin 02), 

(22) 



with ai and 012 defined in Eqs. ([SHZl)- Therefore, 
the coupUng of 02 x X increases as one increases the 
Higgsino components of the LSP. Their coupUngs to 
gauge bosons are 



a2{p)H+{p')W~ 
a^ip)H+{p')W- 
02 (p) hiip') : 



i g2 sin 9a 



2 

i (?2 cos t 



{P " P')t^ : 



1 Cf 

- sin Qa cos (a + /3) (p' - p)^ 



ai(p) /ii(p') : cos6'a cos (a + /3) (p' - p)^ w 



V2 



r 



ig Aff sin 6*^ sin 4/3 , 



V2 



i .g cos Ba sin 4/3 , 



a/2 



(23) 
(24) 
(25) 

(26) 



The coupUngs among 02, a\ and h\ depend on 
the diagonahzation of the CP-even scalar mass ma- 
trix, which may have significant one-loop contribu- 
tions. For simplicity, we use the tree-level results in 
Eq. (fTTj) to obtain analytic formulae. We arrive at 
the following dimensional couplings to the leading 
power in ti/s 

axaxhx : Wa^a^h^ = 0{v^ j s^) , (27) 
02 ai hx ■■ Wa^aihi = -V2KfJ- + C(A^ u^/ s) (28) 

Here the approximation is valid for small values of 
A\ and A^. We also need the main decay channels 
of hi and 02- For the CP-even particle hi with a 
mass below 2Mw, it mainly decays into 2 oi's or 2 
6's with the decay widths calculated as following: 



T{hi 2ai) 
T{hi ^ b + b) 



1 



52nMh, 
3 Af/ii / mi, 

8tt \y/2v 



0(«Vs4), (29) 
2 



(30) 



Thus for u/s ^ 1 the decay hi ^ 2ai is suppressed 
in favor of hi bb. Note this suppression is directly 
connected to the small values of A\ and A^. For 
large values of A\ and A^, the hi 2ai decay 
would be the dominant one p7| . 

Similarly for 02 , if its mass is below twice the top 
quark mass, the leading two decay channels are 



r(a2 hi ai) 
T{a2 ^ b + b) 



K A 
32^ 



Ma, sin 2^, 



3 f "nil, tan /3 sin 6a 



iTT V 



V2v 



(31) 



(32) 



The bosonic decay channel can be dominant for 
modest values for k and A. However, if the 02 mass 
exceeds twice the top quark mass, the decay channel 
into tt opens; 



r(a2 t + t) 



3 ( rat sin B ^ 



Stt VV2i; tan/3 



4m( 

02 



(33) 
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cays to three pions are suppressed by the three-body 
phase space. In the foUowing, we will approximate 
Oi — > 2/x as 100% for case 1. In the second case, 
because the couplings of ai to fermions are propor- 
tional to fermion masses, we anticipate that ai de- 
cays mainly to 2 t's. We do not consider the case 
that ai mainly decays into two electrons because of 
stringent constraints from the beam-dump experi- 
ment at CERN [13. 

To find the interesting parts of the parameter 
space, we use the program NMHDECAY \15^ for nu- 
merical checks. There are many experimental con- 
straints considered in NMHDECAY such as various 
Higgs searches at LEP, 6^57 and T(IS') ai 7. 
The model points presented in this paper pass all 
of the constraints embedded in NMHDECAY. Fur- 
thermore, we will discuss updated constraints on 
T(3S') ^ fli 7 — > fi^ ^~ -f from BaBar and searches 
for dimuon resonances at LEP and Tevatron in sec- 
tion |VlIl 



tan/3 


A 




Ax 








Ml 


M2 


3.1 


0.24 


0.194 


-0.05 




-0.273 


-190 


178.5 


200 




Ma, 


Ma2 


Mn, 




Mh^ 


M„± 


m^± 
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A a2 


161.8 


0.81 


320.1 


114.7 




297.6 


325.4 


175.6 


0.22 


Br(/ii - 


bb) = 78.3% 




Br(;ii ^ 


tt) 


= 8.1% 


Br(/ii ffliai) 


= 1.0% 




Br(ai ^ 


~ 100% 
















Br(a2 ^ 


aihi) = 68.0% 




Br(a2 ^ 


bb) 


= 22.5% 


Br(a2 -» Zhx) 


= 5.4% 




X 


= -0.595 (-iAi) 


+ 0.347 (- 


-iAa) + 0.599 {ipl 


) + 0.404 i^l) 


~ 0.061 {iPs) 


cos 9 J 


1 = 0.12 



TABLE I: ^-favored model point. Masses are in GeV. 

I 

For both cases, in order to isolate the dark matter masses for sleptons, 1 TeV for squarks, 1 TeV for 
discussion, we choose the less relevant soft terms to gluino and —2.5 TeV for all the A-terms in the quark 
be heavy. For example, we choose 500 GeV soft and lepton sectors. We have used the updated top 



and this would become the dominant decay for a2 
assuming k , A < 1 . 

D. Spectrum from numerical calculations 

In this section we find NMSSM model points that 
can provide the neutralino as a DM candidate to ex- 
plain PAMELA. There are two relevant possibilities 
depending on the mass of ai: 

1. The mass of the lightest CP-odd particle ai is 
in the range (2m^, 1 GeV) . 

2. The mass of the lightest CP-odd particle ai is 
in the range (2 rrir, 2 rrib) . 

In the first case oi mainly decays to two muons 
because decays to mesons are kinematically sup- 
pressed. Because ai is CP-odd, decays to two pi- 
ons are forbidden due to CP symmetry, while de- 
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quark mass = 173.1 ± 0.6 ± 1.1 GeV [T^, which 
has a significant correlation to the Higgs boson mass 
for small tan /3 as considered in this paper. 

For the muon- favored case, we choose values for 
other parameters in the NMSSM as in TableHl which 
also shows the relevant spectrum and branching ra- 
tios of the light scalars. The default lower limit on 
the ai mass is 1 GeV in NMHDECAY. One has to 
change the code file named mhiggs . f to obtain an 
ai mass below 1 GeV. As can be seen from Table HI 
if a2 can be produced from w annihilation through 
the resonance effect, the final products of the DM 
annihilation mainly contain ai + hi. The ai decays 
into two muons to provide the positrons needed to 
explain the excess at PAMELA. The dark matter 
mass is mainly controlled by the parameter Mi in 
this model, which is chosen to have the lightest neu- 
tralino mass below the top quark mass. Otherwise, 
02 will decay into t i with a significant branching ra- 
tio, and the neutralino annihilation produces a lim- 
ited amount of positrons and lots of hadrons, the 



hadrons being disfavored by the null antiproton ex- 
cess at PAMELA. It is intriguing that a combination 
of NMSSM and PAMELA results forces us to have 
a dark matter mass below around 170 GeV. 

It is technically natural to have Ma, below 1 GeV 
for tiny values of A\ and as reported here, since 
the U{1)fi symmetry protects its mass (one can also 
use the 11(1) po symmetry to obtain a light pseu- 
doscalar, see [2^ for example). Notice that the 
branching ratio of hi — > aioi is below 1.2% to sat- 
isfy the current null results of searches of oi in the 
channel hi — > aiai ^ 4 /i at DO (see Section IVlIB[) . 

For the tau-favored case, we list the values of 
model parameters, spectrum of particles and inter- 
esting branching ratios of light scalars in Table [ill 
The current direct searches only impose mild con- 
straints on the model parameters. Therefore, we 
choose one representative point in the parameter 
space to have hi aiai At as the main decay 
channel of hi, and hence to have six r's in the final 
state of dark matter annihilations. 
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Mh, 
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Br(/ii - 


aiai) = 92.2% 
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^ hb) 


= 7.6% 










Br(ai - 


+ T+T-) = 85.0% 
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99) 


= 7.7% 






Br(ai cc) 


= 5.4% 


Br(a2 - 


aihi) = 95.5% 




Br(a2 


bb) 


= 2.8% 






Br(a2 — » 01^2 


) = 0.7% 




X = 0.968 (-lAi) 


+ 0.003 ( 


-1X2) - 0.236 (iP'i) 


- 0.080 M) 


+ 0.043 {iPs) 







TABLE IL r-favored model point. The masses are in 
GeV. 



III. DARK MATTER ANNIHILATION ter annihilation products. This goal can be achieved 

CROSS SECTION 



Our first goal is to look for parameter space in 
the NMSSM having leptons as the main dark mat- 
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if the lightest CP-odd scalar is an intermediate an- 
nihilation product and subsequently decays to two 
taus or two muons from kinematic constraints. The 
other goal is to have a large dark matter annihilation 
cross section to explain the size of the PAMELA ex- 
cess. The annihilation cross section can be enhanced 
through s-channel resonance effects. Assuming the 
lightest neutralino is the dark matter candidate in 
the NMSSM, the heavier CP-odd scalar 02 is the 
only particle which can play this role. The CP- 
even scalars are ruled out by the CP symmetry, be- 
cause the initial state with two identical Majorana 
fermions is CP-odd. The lighter CP-odd scalar ai 
is far below the necessary mass region for the reso- 
nance effect. 

Having established XX ~^ 0,2 X as the main 
dark matter annihilation channel, the main annihi- 
lation products are equivalent to the decay products 
of a2. The main dark matter annihilation products 
are shown in Fig. [1] for the two cases exemplified by 
Table Hand Table HIl 




FIG. 1; The Feynman diagram of the dominant neu- 
tralino annihilation channel. The upper panel is for the 
muon-favored case, while the lower panel is for the tau- 
favored case. 

To calculate the dark matter annihilation cross 
section, we will concentrate on the main annihilation 
channel XX ^ '^2 ~^ 0,1 hi, which is true for both 



two cases. The annihilation rate is 
cr{xx ^ 02 ^ ai hi)vdm 

^ 9I2XX '^aia2hi _ Ml \ 

len [{Ml ~ AmlY + M^] \ Ami)' 

(34) 

up to (^{v^^/c^) where the average dark matter 
speed is Vdm/c ~ 10""^ in the galactic halo. Here 
the mass of oi is neglected and the coupling Wa^a^hi 
has mass dimension one. The co-annihilation effects 
can be neglected here, because other superpartner 
masses are much larger than the LSP mass. For 
the parameter points in Table |T] and Table [iTl we 
have <j{xX ~^ aihi)vdm approximately 135 pb-c 
and 301 pb-c, respectively. Also from Table |T] and 
Table [III the width Fqj of 02 is below 1 GeV, so 
the width part in the denominator of Eq. (j34p can 
be neglected for \Ma2 — 2m^| > 1 GeV consid- 
ered here. Therefore, a large dark matter annihila- 
tion cross section can easily be obtained for the case 
of < mt- On the contrary, if > mt, the 
decay channel of a2 — > tt is open and the decay 
width of 02 is of order 10 GeV. Then the dark mat- 
ter annihilation cross section is limited by the width 
part in Eq. p4)) . generically below 10 pb-c and not 
large enough to explain the PAMELA data. Thus 
we have two reasons to believe that the lightest neu- 
tralino mass is below the top quark mass: one is 
to have dominantly leptonic annihilation final states 
and the other one is to have a large annihilation cross 
section. 

Our annihilation rates are by far larger than the 
necessary one (~ 1 pb-c) to satisfy the dark matter 
thermal relic density. One possible explanation for 
this discrepancy is that the dark matter is nonther- 
mal. For example, other long-lived particles like the 
gravitino or moduli can decay into the LSP at a later 
time [ld|. As argued in ^1U\, a long-lived modulus 
field with a lifetime shorter than one second can nat- 
urally appear in the anomaly-mediated SUSY break- 
ing model. The main dark matter relic abundance 
will be determined by the moduli annihilation cross 
section, which in principle can be smaller than the 
one for the neutralino and provide the observed dark 
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matter relic density at the current time. 



IV. POSITRON EXCESSES FROM 
NEUTRALINO ANNIHILATION 

We now calculate the positron spectrum out 
of dark matter annihilation and compare it with 
the PAMELA results to determine the PAMELA- 
favored parameter space in the NMSSM. We 
first discuss the source term of primary positrons 
from DM annihilations, propagation of cosmic ray 
positrons and then the positron fluxes measured at 
PAMELA. We also consider the charge-dependent 
solar modulation on the positron fraction spectrum. 



A. The source term for primary positrons 

For the two cases we are considering in this pa- 
per, the source term for primary positrons can be 
generally written as: 

,(x,i?) = j (35) 

Here the overall factor 1/2 is from the Majorana 
property of the neutralino DM candidate; (cr v) is the 

I 



Here V{Ei — > Ej) denotes the probability of a parti- 
cle i with energy Ei decaying into a particle j with 
energy Ej; the factor of 2 in the expression is be- 
cause hi decays into 2 ai's. We neglect the positrons 
and electrons from the b quark decays although the 
hi mainly decays into two b quarks in the muon- 



thermally averaged annihilation cross section and to 
a good approximation can be replaced by the for- 
mula in Eq. ([M]) ; dNg^+ / dE^+ is the energy spec- 
trum of positrons; p(x) is the dark matter distri- 
bution inside the Milky Way halo. For the dark 
matter distribution we use either the Navarro, Frenk 
and White (NEW) profile [IH or the cored isother- 
mal (ISO) profile |22]; the use of other profiles may 
change discussions in this paper, especially for the 
gamma ray spectrum. The NEW and ISO profiles 
are parametrized as: 



PNFw(n = Pq [ — \ 

\ r / \ rs + r J 

Piso(r) = Pe [^^) ' (36) 

with r, = 20 kpc (NEW), 5 kpc (ISO) is the 
radius of the central core; r© = 8.5 kpc is the 
galactocentric distance of the solar system; pq = 
0.3 GeV cm~'^ is the solar neighborhood DM den- 
sity. 

The positron energy spectrum function for the n 
case can be generally expressed as 



(37) 

I 

favored case. Because positrons and electrons from 
b decays are relatively soft, including them makes 
only a slight change in the spectrum below 10 GeV, 
where the background is anyway dominant. 

Erom two x's annihilating into ai and hi, and 



dE, 



■ViEa,^ E.+ )V{E^+ ^ E,-,) 



2 J dEh, dEa, dE^+ ^ V{Eh, ^ EajViEa, -> E^^)V{E^+ ^ E, 
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neglecting the mass of ai, we have 



dEa, 

dNh, 
dEh, 



Br(a2 hiai)6 Ea^ - (m^ 



Br(a2 hiai)S Eh^ - (m^ 



Ml 

4 TTlv 



Since the emission of muons in the ai rest frame and 
the emission of oi's in the hi rest frame are isotropic, 
we have 

ViEa.^E^^) « -^n{Ea, -E^+), 



V{E^,^Ea,: 



Br(/ii 



ai ai) 



E, 



hi 



(39) 



with 7i(a;) is the heavy-side function and the masses 
of fi^ and ai are neglected in the approximation 
formula. Neglecting the muon polarization, the 
positron energy probability in muon decay has the 
following analytic form [3] 



Ee+) 



1 



3E^+ 



[5 - 9x^ + 4x^] 



xh{E.+ - E,+ ), 



(40) 



with X = E^+/E^+. Here the functions V{Ei — > 
Ej) are normalized such that J dEjV{Ei ^ Ej) — 
Br(*^j). 

The r-favored case is more complicated than the 
/x-case, because other than leptonic channels r can 
also decay to various charged mesons, which decay 
eventually to electrons and positrons. In order to ob- 
tain the electron energy dependent probability from 
r decays, we use PYTHIA [23|, which calls the pro- 
gram TAUOLA [24|, to simulate the inclusive elec- 
tron/positron energy spectrum from both direct and 
indirect r decays. As can be seen from Fig. [51 once 
the energy of the (unpolarized) r is much larger than 
its mass, the fractional electron and positron en- 
ergy spectra are independent of the energy of the 
r. The large fluctuations for £'g-/£'^+ > 0.5 in 
the lower panel of Fig. [5] are due to limited Monte 
Carlo statistics. The tiny fraction of e~ out of r"*" 



0.1 ' 



0.01 



-a 



0.001 



(38) 
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FIG. 2: Upper panel: the positron energy spectrum from 
inclusive tau decays. The red and the black points are for 
10 GeV and 300 GeV r+ energy, respectively. The blue 
line is from a fitted analytic function in the text. Lower 
panel: similar as the upper panel, but for electrons. On 
average, 1.16 positrons and 0.16 electrons are generated 
from one decay. 



for E^- /E^Jr > 0.5 can be understood from the fact 
that only high multiplicity final states can contain 
an e~ with a different charge from r^. 



For these e"*" and e spectra from inclusive r"*" 
decays, the following fitted functions provide good 
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agreement: 



1 



-97.716 +223.389 2;'' 



X e 



Er+ 

193.748 2;^ +82. 595 - 22. 942 a; + 2.783 ^j-^ _ 
1 



E^ 



-15.5752;'' + 15.792;^ 



18.083 2; + 0.951 



E^+ 
nil - a;). 



(41) 



with X = E^,± /Et+. For each r+, on average, there 
are j dE^+V{E^+ E^+ ) m 1.16 positrons and 
/ dE^-V{E^+ Ef.-) « 0.16 electrons produced 
in the decay. When ai decays into t+ + r^, there 
are on average 1.32 positrons out of a single ai de- 
cay. Therefore, we insert Eq. (|4T|l into the analogous 
formula in Eq. (|37p for the r case. 



B. Propagation of positrons and electrons 

The propagation of positrons and electrons in the 
galactic medium is described by the following trans- 
port equation: 



dN 
'dt 



V • [K{yL,E)\/N] 



d_ 

dE 



[b{E)N]=q{^,E){^2) 



where A^(x, E) denotes the positron number den- 
sity per unit energy; qi^, E) is the positron source 
term; K{x,E) = Kq{E / EqY is the diffusion con- 
stant with Eq = 1 GeV; b{E) = E^/{EoTe) is 
the positron energy synchrotron and inverse Comp- 
ton loss rate with te = 10^^ s. The diffusive halo 
is modeled as a cylinder with radius — 20 kpc 
and the vertical direction z inside {—L, L). The half 
thickness is not constrained and varies from 1 to 
15 kpc. We will consider three different parameter 
points for the cosmic ray propagation model in Ta- 
ble mil Those sets of propagation parameters are 
compatible with the secondary /primary test for the 
secondary and primary antiprotons [26'|. Assuming 
a time-independent state and considering the diffu- 
sion constant as space independent, we have 



AN 



d_ 
'dE 



E' 



Eqte 



-N 



g(x,i?)(.43) 



Model 


S 


Ko [kpcVMyr] 


L [kpc] 


M2 


0.55 


0.00595 


1 


MED 


0.70 


0.0112 
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Ml 


0.46 


0.0765 
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TABLE III: Three combinations of cosmic ray propaga- 
tion parameters, which give the minimum, median and 
maximum positron fluxes [25j . 



Defining a pseudo-time and a relative pseudo-time 
between the source point and the observation point, 
respectively, as 



t{E) = TE J— — ^ , 

f{E,Es) = i{E) - i{Es), (44) 



the characteristic diffusion length in the radial di- 
rection is 



\d{E,Es) = y/AKofiE,Es). 



(45) 



In Ref. [25j, an analytic solution for the positron 
flux on the Earth has been obtained and has the 
following form: 



€AE) = ^ r dE^q{rQ,E,) 
47r Je 

xr]{XD{E,Es)) . 



teEq 
E^ 



(46) 



Using the Bessel expansion method, the halo inte- 
gral rj, which is the volume integration of the Green 
function of Eq. has both radial and vertical 

expansions: 
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OO CXj 

Vi^o) = X! X! M^irQ/rs) sin( — )exp 

i—l n—1 

with 

-Ri.n = -TT — 7^-^ / drrJo[x„r/rs)— / dzsin(-— ) . (48) 

Jiix^yrj Jo L 2i \^ P0 J 




Here Jfe is the Bessel function of the first kind and 
Xi is the i-th root of the Jq Bessel function. In prac- 
tice, one can use the numerically fitted functions in 
Ref. [131 to speed up the numerical calculations. For 
example, we use the following numerical function for 
the NFW dark matter profile and the M2 propaga- 
tion model: 



Tj{Xjj) = 0.5 + 0.774 tanh 
X -0.448 exp 



0.096 



0.211 

{i - 192.8)2 
33.88 



0.649 



(49) 



with £ = logio(A£)/kpc). 



C. Solar modulation 

PAMELA has measured the positron over electron 
fraction with energy below 10 GeV, and obtained a 
spectrum significantly below the background fitted 
from other cosmic ray experiments. One possible 
explanation of this discrepancy is due to the charge 
sign dependence of the solar modulation. The mag- 
netic field of the solar wind is dominated by the 
dipole term, and the projection of this dipole on the 
solar rotation axis can be either positive or nega- 
tive, called and A~ states, respectively. At each 
sunspot maximum, the dipole reverses its direction 
and leads to a periodic function for the dipole mag- 



netic field with a roughly 12 year period. ^ 

Using two functions c+{E) and C-{E) to model 
the solar modulation, we have the observed positron 
fraction on the Earth as: 



c±{E) $,+ (£;) 



c±{E)^,+ {E) + c^{E)<P,-{E) 



(50) 



with -I- for the solar system in the A'^ state and 
— for the A^ state. Here $e± (E) denote the 
positron/electron fluxes. One notices that only the 
ratio of c+{E)/c-{E) is relevant for the positron 
fraction. The ratio of the total electron flux in the 
A~^ cycle to total electron flux at a similar phase in 
the A^ cycle is 



RiE) = 



c+{E)<i>,+ {E) + c^{E)<i>,-{E) 
c-{E)<P,+ {E) + c+iE)^,-{E) 



(51) 



This ratio can be fitted in principle from the data 
shown in Fig. 5 of Ref. [28|. However, since the 
lower energy bins of PAMELA (where the contribu- 
tion from dark matter is negligible) are measured so 
precisely, one can use the data from PAMELA {A~ 
cycle) and the known background [1^ (without an 
assumption of solar modulation) to determine this 
ratio. The following fitted formula can achieve the 
goal: 

R{E) = max [min (0.48 + 0.2 log(£;/GeV), 1.0) , 0.2] (^2) 



1 Figure 1 of Ref. ^ 
in the electron flux. 



shows clear evidence for this behavior 
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which wiU be used in our following analysis about 
the comparison of model predictions and PAMELA 
data. The general form of this expression is sug- 
gested by Ref . js^ . Eventually one needs to under- 
stand or to calculate this ratio function or c±{E) 
from first principles. Due to the magnitude of the 
solar magnetic field, only positrons/electrons with 
energy below about 10 GeV can be influenced by 
the solar wind. R{E) is a monotonically increasing 
function and saturates at unity above 7 GeV. 

In terms of R{E) and F{E) = 
$e+(-E)/(<i>e- (i?) -t- <^p+{E)), which is the positron 
fraction without solar modulation effects, one has: 



PAMELA results, we first discuss the electron and 
positron backgrounds from standard astrophysicial 
processes. Background positrons are mostly secon- 
daries originating from spallation processes of cosmic 
rays, mostly primary protons, off the interstellar gas, 
thus mainly occurring in the galactic disk. The pri- 
mary electrons are mainly produced by shock wave 
acceleration in supernovae. To simplify our com- 
parison of model predictions to experimental data 
(without the solar modulation effects), we use the 
following numerically fitted formulae for the back- 
ground 1301 ■ which agrees with the full results calcu- 
lated in I29II: 



F^{R + 1) - FR 
2F - 1 



= FHR + i)-F 

R{2F - 1) 



(53) 



D. PAMELA positron excess 

Before we compare the positrons from the 
neutralino annihilations in the NMSSM to the 



,bks 



^sec,bkg 
^scc.bkff 



0.16 £; 



-1.1 



1 + llE^' 



i- 3.2£;2 i5 
0.70 



580£;4-2 ' 

_x T-_650J2.3 + i5ooi;4.2 ' (5"^) 
in GeV cm s sr and with E in GeV. Com- 
bining the background and signal positrons, we have 
a general formula for the positron fraction 



1 + nOE^-^ + 600E^-^ 
4.5 E" '^ 

1 



F^ = 



sec,bkg 



c+. 



prim,bkg 



scc,bkg 



+ 0lig, 

is) + 



(55) 



scc.bkg 



assuming that the PAMELA data were taken when 
the sun is in the cycle. 

Using the model point and the masses of x, hi and 
ai reported in Table [H we have the positron excess 
for the /^.-favored model point in the NMSSM shown 
in Fig. El To generate the plot in Fig.jSj we have used 
the M2 propagation model from Table IIIIl which 
provides a best fit to the PAMELA data. The other 
two propagation models generate a flatter curve than 
the M2 model. This is because as the thickness of the 
diffusive halo decreases, the positrons detected on 
the Earth originate from a nearby region (the char- 



acteristic propagation distance Ad decreases), and 
hence low-energy positrons are less likely to reach 
the Earth. This leads to a steeper spectrum for the 
M2 propagation model. Since all three propagation 
models are supported by the N-body simulation, the 
combination of the NMSSM and the PAMELA data 
(assuming the NMSSM interpretation was confirmed 
by, e.g. collider discoveries) could help to determine 
a correct galaxy model. 

In Fig. [3] we note that the NMSSM model is 3 
sigma below the PAMELA data point in the highest 
energy bin. However the overall agreement is quite 
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good. To estimate the goodness of fit conservatively, 
we ignore the PAMELA data below 7.4 GeV, which 
should have negligible contribution from the dark 
matter annihilation, and calculate the averaged 
for the 8 bins above 7.4 GeV: 




1 2 5 10 20 50 100 



Energy (GeV) 

FIG. 3: The positron fraction from the /x-favored model 
point in the NMSSM. The sohd black hne is the back- 
ground without considering solar modulation. The dot- 
dashed black line is the background with the solar mod- 
ulation effects. The dashed blue line is the positron 
fraction from neutralino (m^ = 162 GeV) annihilations 
plus the modulated background. The red points are data 
from PAMELA with one standard deviation errors. The 
dark matter annihilation cross section is 6.0 x 10~^* cm^ 
s~^. The M2 propagation model is used here. 

If the neutralino annihilation in the NMSSVI is the 
explanation for the PAMELA data, the rising fea- 
ture of the positron fraction spectrum should end at 
around 70 GeV. This is a dramatic prediction for fu- 
ture PAMELA results. This predicted turnover fol- 
lows uniquely from the requirement in the NMSSM 
model that the dark matter LSP mass cannot exceed 
the top quark mass. 

For completeness, we also show the electron plus 
positron energy spectrum from dark matter annihi- 
lation in Fig. m We also include the latest results 
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FIG. 4: The electron plus positron spectrum from the 
/u-favored model. The solid black line is the background. 
The dashed blue line is from the neutralino annihilation 
plus the background. The red crossed points are the 
results from Fermi Large Area Telescope (Fermi LAT) 
with the gray band for systematic errors. The dark mat- 
ter annihilation cross section is 6.0 x 10^^'' cm^ s~^. The 
M2 propagation model is used here. 

from Fermi LAT in the red and crossed points [3l| in 
Fig. m There is an additional uncertainty from the 
LAT energy scale, which can shift the whole gray 
band by 5% (up) to 10% (down) and is not shown in 
this figure. In this figure the electron and positron 
fluxes from the background in Eq. ([5^ plus signal 
have been normalized to agree with the first bin of 
Fermi LAT. The positron fraction predictions from 
PAMELA are unchanged by this manipulation. 

The agreement between the predicted elec- 
tron-|-positron spectrum and the Fermi LAT data 
is about the same with and without adding the 
NMSSM signal. The generally poor agreement 
should thus be attributed to a defect in our un- 
derstanding of the cosmic electron/positron back- 
ground. More generally, we conclude that a neu- 
tralino annihilation explanation for PAMELA is con- 
sistent with the Fermi LAT results, as long as the ex- 
tra contribution to the electron-f-positron spectrum 
is within the Fermi LAT errors. 

Having discussed the /j,-favored point in the 
NMSSM, we also report the results for the r-favored 
point. It turns out that although r-favored points 
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can be found easily in the parameter space of the 
NMSSM, they provide a worse fit to the PAMELA 
data. This is mainly because the positrons from t+ 
decays are softer than the positrons from /i+ de- 
cays. This fact can be seen from the lower panel of 
Fig. [5l where a comparison between /x and r cases 
is shown. Only taking the PAMELA data above 
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V. PAMELA ANTIPROTON 

Since there is no excess of the antiproton frac- 
tion observed at PAMELA, this imposes constraints 
on the antiproton production cross section from 
dark matter annihilation. Specifically to the muon- 
favored case in the NMSSM, the dominant source of 
hadronic production is from ft-i — > 6 6 and 02 — > 66. 
Similarly to the calculations for the positron frac- 
tion spectrum, we first use an analytic function to 
fit the fragmentation function of 6 6 to antiprotons. 
The antiproton energy spectra for two different 6 6 
center of mass energies are extracted using PYTHIA 
and shown in Fig. [HI We use the following fitted 
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FIG. 6: Tlie 6 6 to antiproton fragmentation function. 
Tlie red and black points are for 20 GeV and 300 GeV 
66 center of mass energies. The blue line is from a fitted 
analytic function described in the text. The fragmenta- 
tion function is extracted using PYTHIA. 



FIG. 5; Upper panel: the positron fraction from one 
r-favored model point in the NMSSM. The dark matter 
{m-^ = 160 GeV) annihilation cross section is 9.0 x 10~^* 
cm'^ s~^. The M2 propagation model is used here. Lower 
panel: a comparison of positron energy spectra from r"^ 
decays and from n'^ decays. 



7.4 GeV into account, we calculate the average 
for the fit of this r-favored model point to PAMELA 
data as ~ 3.4. Because of this large x^, we wiU 
concentrate on the /i-favored model from here on. 



function later to replace the numerically simulated 
antiproton energy spectrum from hi decays: 



X e 



bb Tp) ^ -^+525.86. 
-331.51 x^+96.71 -28.56 a;+2.755 



(57) 



where x — 2 {Tp + mp)/Eh-^. As is customary in 
antiproton cosmic ray analyses, we substitute the 
proton energy by its kinetic energy Tp = Ep — rrip 
in the following. The antiproton spectrum from 02 
has a same formula replacing E^-^ with Ea2 ■ 
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The source term for antiprotons from DM an- 
nihilations is similar to Eq. (|35p but substituting 
dN^+/dE^+ with 



(riH + 4^/^ riHc) Cp?" vp with nn ~ 1 cm 



dTp 



— — Br(/ii 



hh 
bb 



Tp) 



dE,. 



bb) 

dNg 
dEa 



Br(a2 



> bb) 
(58) 



However, due to the fact that 3> rrie, the en- 
ergy loss term for antiprotons can be neglected. The 
steady diffusion equation for antiprotons is [26l |: 

-Kp{T)AN + |-(sign(z)KiV) 



dz 
2hS{z) 



,N = 9(x,T)(.59) 



Here iV is the number density of antiprotons per 
unit energy and Kp{T) = KoP{p/GeY)^ with /3 
and p are the antiproton velocity and momentum. 
The second term is related to the convective wind, 
which has a direction outward from the galactic 
plane and represents the movement of the medium 
responsible for the antiproton diffusion. The ve- 
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TABLE IV: Three combinations of cosmic ray propaga- 
tion parameters, which give the minimum, median and 
maximum signal antiproton fluxes. 



locity Vc is assumed to be constant and has val- 
ues shown in Table IIVI for different propagation 
models. The third term represents annihilations of 
antiprotons and interstellar protons in the galac- 
tic plane with a thickness ft, = 0.1 kpc. The an- 
nihilation rate between antiproton and protons is 



The cTp^n as a function 



and riHo ~ 0.07 cm 

of antiproton kinetic energy is given in |32| and 
dl]: 661(1 + 0.0115 T-O-^^" - 0.984T° °^^^) mb for 
T < 15.5 GeV and 36T'°-^ mb for T > 15.5 GeV. 
This annihilation process is dominant at low energy 
and leads to a decreased flux of antiprotons with 
low energy. In our analysis, other non-annihilation 
interactions between antiprotons and the interstellar 
medium are neglected. These effects are not impor- 
tant for the antiproton flux with energy above a few 
GeV. 

Similar to the positron case, the diffusion equation 
for the antiproton can be solved analytically and has 
the following concise form as its solution [2a |: 

For the M2 propagation model and the NFW dark 
matter profile, we numerically fit an analytical for- 
mula for R{Tp) 

i?(T)/Myr = IQI-^SS + 0.0542 logio(T) -0.265 logio(T)=* 
^ -^Q + O.OSQT logio(T)^ -,0.00575 logio(T)'' 

(61) 

with T in GeV. To obtain the antiproton flux ob- 
served on the Earth, we need to take the solar mod- 
ulation effect into account. The solar modulation 
effect is represented by a parameter 0, which is 
500 MV for minimum solar activity when PAMELA 
was taking data. The energy spectrum of antipro- 
tons on the Earth is 



'PfiTp) = (b^iTp 



{Tp + nip)^ - ml 



(Tp 



(62) 



with Z = 1 the electric charge of the antiproton. 

The background for primary protons can be ex- 
trapolated from other cosmic ray experiments. Us- 
ing the data from AMS [13] , we arrive at the follow- 
ing fitted analytic function to describe the primary 
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proton background: 

0BG,AMS ^ji^ _ g-0. 00097 log'' T + 0.012 log* T 
^ g + 0.014 log-' T- 0.382 log^ T - 0.828 log T + e.SSj-^g-j 

with T in GcV and the flux in m^^ s^^ sr~^ GeV"^. 
Similarly, from CAPRICE98 [35] data, we have 

^BG, CAPRICE j-jn^ _ g-0.0005 log^T + 0.005 log'' T 
^ g+0.019 log'^T- 0.433 log^T- 0.882 log T + 6. 89 ^g^-j 

The primary proton fluxes measured at AMS and 
CAPRICE can deviate from each other by a differ- 
ence as large as 20%. We include both results in 
our analysis to encode uncertainties of our primary 
proton flux background. The secondary antiproton 
background can be found in the detailed analysis 
in [201 and fitted in [131 as: 

^BG^jn^ ^ qO. 028 log*(,T- 0.02 log^^T 

^ 10 " ^°<^"> ^ + ^-""^ '°Sio T - 1.64 ^ (-gg^ 

A comparison of the antiproton/proton flux ratio 
from the neutralino annihilation and the PAMELA 
data is shown in Fig. [T] Noticing that the M2 
propagation model provides a better fit to the 
PAMELA positron data, the model point in Ta- 
ble HI is marginally allowed by the PAMELA an- 
tiproton data. To quantify the discrepancy between 
the model prediction and the PAMELA observed 
data, we calculate the by including all bins. We 
also calculate the between the background and 
PAMELA data. Using the background extrapolated 
from AMS, we have to be 1.6 for the DM predic- 
tion plus the background and 3.6 for the background 
only. While, using the background from CAPRICE, 
we have x^ to be 4.0 for the DM prediction plus the 
background and 2.2 for the background only. We 
conclude that the model point in Table [T] is allowed 
by the PAMELA antiproton data, taking into ac- 
count the large uncertainties of the background pri- 
mary proton fiux. 

Although the branching ratio Br(a2 6 6) is 
smaller than Br(a2 — s- hi ai)Br(/ii 6 6), the center 
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FIG. 7: The antiproton-to-proton flux ratio as a func- 
tion of kinetic energy. The black solid and dashed lines 
are the backgrounds from AMS and CAPRICE, respec- 
tively. The thick green solid line is from the dark mat- 
ter annihilation plus the background (AMS) using the 
M2 propagation model. The thick blue dashed line 
is using the background (CAPRICE). The red crossed 
points are data from PAMELA. The neutalino mass is 
— 162 GeV and its annihilation cross section is 
6.0 X 10-2'* cm^ s-\ 



of mass energy of 66 directly out of 02 is approxi- 
mately twice of the center of mass energy out oi hi. 
Therefore, the large antiproton fraction for the ki- 
netic energy above 10 GeV mainly comes from dark 
matter annihilating directly into 6 b. To suppress the 
antiproton flux more efficiently, one could change the 
model parameters to suppress the branching ratio of 
02 to 66. One simple way to do this is to reduce 
tan/3. However, by doing so, the branching ratio of 
hi ai ai is also increased, creating a tension with 
the upper limit from DO described in Section fVlI Bl 



VI. GAMMA RAY FLUXES FOR FERMI 
LAT 

If the ^-favored model points in the NMSSM are 
the correct explanation of the PAMELA positron 
excess, there will be lots of associated gamma r ays 
generated. Existing gamma ray data from HESS [37[ 
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and EGRET [38[ can in principle impose constraints 
on the dark matter annihilation cross section to elec- 
trons and positrons. Also, the recent and upcom- 
ing gamma-ray flux data from Fermi LAT have a 
smaller statistic uncertainty and can be used to test 
the PAMELA-favored NMSSM model. 

The differential gamma-ray flux from the dark 
matter annihilation has the following general for- 
mula: 



dndEl 



1 



{(TV 



2 Aw ml dE-y 



p^{r)dl{il)). (66) 



Here — P + — 2lrQ cos ip with ip as the an- 
gle between the line of sight and the galactic plane. 
One can separate the astrophysical uncertainties by 
introducing the quantity 



J(V) 



1 



p\r)dm. 



(67) 



Performing the solid angle integration, the differen- 
tial gamma-ray flux is 



d<^^ 
dK 



1 r-Q Pq (av) dN. 

2 An ml 'dE. 



An , (68) 



with J{An) = (l/Arj) J^^J{i;)dn and An = 
27r(l — cosip) for the region around the galactic 
center. For example, we have J(Ail) An w 1 for 
An = 10-3 and J{An) An « 0.1 for An = 10^5 sr 
using the NFW dark matter profile [soj . 

The gamma-rays from the /x-favored model point 
of the NMSSM have two sources: one is related to 
the muons in the final state and the other one is 
related to the bottom quarks in the annihilation final 
state. Altogether, we have 



dN^ 



Br(a2 ai hi) 



Br(a2 ^ ai hi) 



dNi;{Ea,) 

dE^ 
dN^jEh,) 

dE^ 



Br(a2 bb) 



dN^jEg,) 
dE^ 



(69) 



There are two processes to generate gamma-rays as- 
sociated with the muon final state. One is through 



final state radiation (FSR) and the other one is from 
the radiative muon decays into photons. For the fi- 
nal state radiation, we have 



dE., 



1 + (1 - x)^ 



2 

Eai 

log 



dx 

X TT 



Mlil-x) 



(70) 



For the radiative muon decays: ^ — s- e v^Vel and 



e+ VeV^^, one has [40 

dNt;^„An{Ea,) ^ 
dE^ 



E„, 



dx — 

x 



dy-T{y) (71) 



with (3 = y^l — Am'^/m'^ and J-{y) as the photon 
spectrum in the muon rest frame and given by 



- (3-2y + 4y2-2y3)log^ 

d TT y \ m^ 

17 23 y 101 y2 55 y3 
" y ^ ~6 12" ^ 12 

+ (3 -2y + 42/2 _22/3)log(l ~ y)) ,(72) 

with y — 2Ery/mp,. 

The gamma-ray fragmentation function from bot- 
tom quarks in the dark matter annihilation final 
state is simulated using PYTHIA and fitted using the 
following analytic function: 



dN^JE) 



dE^ 



^7.59 - 



3.25 a: -I- 412.5 a:^ - 1297. Tx'' 



E 



xe 



+ 1969.7 2:* - 1137.8 



(73) 



with X = 2E^/E and E as the center of mass energy 
of the b b system. 

Summing up all contributions to the gamma-rays, 
we compare the model predictions with the back- 
ground for gamma-ray energy above 1 GeV, which 
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is fitted by a power-law in Ref. [4l| as 



VlGeV 



2.72 



(74) 



The Fermi LAT collabora- 
tion has already shown a preliminary result for the 
gamma rays from 0° < / < 360° and 10° < |6| < 
20° [ill. Here, I and h are the heliocentric galactic 
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FIG. 8: The gamma ray flux predicted from the neu- 
tralino annihilations in the NMSSM. The neutralino 
mass is = 162 GeV and its annihilation cross sec- 
tion is 6.0 X 10"^'' cm'^ s~^. The solid green line is the 
dark matter prediction using the NFW proflle plus the 
traditional background (shown as the thin black line). 
The dashed orange line is the dark matter prediction 
using the cored isothermal dark matter proflle. The dot- 
dashed green line is the dark matter prediction without 
including the contributions from 6 quarks. 

coordinates. This region of angles corresponds to 
A57 = 0.567 sr. It is easy to calculate JAfi to be 
13.4 for the NFW profile and 5.7 for the ISO pro- 
file. In Fig. [HI we show the observed gamma ray 
fluxes from Fermi LAT together with the data from 
EGRET. As one can see, there is a disagreement 
between those two experiments for photon energy 
above 1 GeV. The predictions from dark matter an- 
nihilation using the model point in Table [J tend to 
agree with the EGRET result if one uses the ISO 
dark matter profile. 



A larger discrepancy occurs for the dark mat- 
ter prediction using the NFW profile. This might 
be reconciled by astrophysical uncertainties like the 
smoothness of the dark matter distribution, which 
affects the necessary dark matter annihilation cross 
section by a factor of a few. Notice also in Fig.[8]the 
dot-dashed green line showing the signal gamma-ray 
contribution from the muon final states only; this 
indicates that the discrepancy is coming from the h 
quarks in the final state. 

We can also consider gamma rays coming from 
the galactic center. In Fig. [9l we compare the dif- 
ferential gamma-ray flux predicted from the neu- 
tralino annihilations in the NMSSM to the flux from 
the background, after fitting the positron fraction 
spectrum of PAMELA. We present the comparison 
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FIG. 9; The photon flux predicted from the neutralino 
annihilations in the NMSSM after fitting the positron 
fraction spectrum of PAMELA. The neutralino mass is 
= 162 GeV and its annihilation cross section is 6.0 x 
10"^* cm^ s"^ An = 2.4 x 10~* (0.5° x 0.5° about the 
galactic center) and JA^ ^ 0.69 (NFW) and JA^i f» 
0.003 (ISO) are used here. 

of the model predictions and the background for a 
0.5° X 0.5° region about the galactic center. Since the 
Fermi LAT will have an angular resolution of around 
0.1° X 0.1° (10~^ sr), the error of their measurements 
in principle is small enough to test the dark matter 
annihilation scenario. However, the model predic- 
tions for the gamma-ray fiux in the galactic center 
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are subject to large uncertainties from the dark mat- 
ter profile, which can bring a factor of few hundred 
difference as seen from Fig. [HI 

Due to kinematical reasons, our dark matter neu- 
tralino mass is less than the top quark mass. There- 
fore, the gamma rays predicted in this model have 
an energy cutoff below the top quark mass. Hence, 
this model automatically evades the constraints from 
HESS, which measured the gamma-ray flux with en- 
ergy above 200 GeV. For EGRET with Aft = 10"^ 
around the galactic center and with the energy range 
of 1 GeV < < 30 GeV, the total background flux 
is 3 X 10~^ cm^^ s~^. The dark matter contribution 
in this model is around 2 x 10~^ cm~^ s~^ for the 
NEW profile and 1.5 x lO^^ cm^^ g-i for the ISO 
profile. This indicates a tension between the model 
prediction with the NEW profile and the EGRET 
data. Again this result is subject to astrophysical 
uncertainties; for a solid angle around the galactic 
center, changing the dark matter profile can intro- 
duce an uncertainty of two orders of magnitude in 
the gamma-rays flux predictions. 



VII. DIRECT CONSTRAINTS 

There are many direct constraints on the NMSSM 
parameter space from LEP, Tevatron, CLEG, B- 
factories and the magnetic moment of the muon. 
Our analysis shows that the most stringent con- 
straints are from exclusive T(3S') decays at BaBar 
and light pseudoscalar searches at DO. Compared 
to the MSSM, the constraints from LEP and other 
searches are less severe, due to non-standard decays 
of the MSSM-like Higgs boson and a new suppres- 
sion factor cos 9a beyond the MSSM. 



A. LEP constraints 

The heavier bosons in the NMSSM, typically 
above 250 GeV, were unaccessible at LEP. There- 
fore, we only consider the two lightest scalar parti- 
cles ai and hi at LEP. There are three main pro- 
duction mechanisms for those two neutral Higgs 



bosons. One is through the Higgsstrahlung process 

e~ — !■ hi Z] another one is the pair production 
process — > ai; the third one is through the 
radiation off a massive fermion: e+ ^ bbai. 

For the Higgsstrahlung process e+ hi Z and 

for an hi mass within the LEP reach, the hi mainly 
decays to 6 6 and 2 oi for the /z-favored and r-favored 
cases, respectively. For the r-favored model point, 
there are many different final states like 4t, 2t 2g 
and so on. Although the LEP bounds on each of 
those channels are weak, the decay mode indepen- 
dent limits impose a bound on the mass of hi as 
> 82 GeV [43*], which is satisfied in the model 
points we have considered in this paper. For the 
/.t-favored case, the lower bound on the hi mass is 
around 114 GeV, which is also satisfied for the model 
point reported in Tabic [H 

The cross section of the pair production process 
e"*" e~ hiQi is proportional to cos^ 9a M^/M^^ 
from Eq. ([26|) . and hence is tiny for cos 0a < 0.3 
and Afa2 > 300 GeV. For = 200 GeV, the cross 
section is calculated to be 2 x 10"'^ fb for the model 
parameters in Table [J Considering the integrated 
luminosity of LEP is below 1 fb~^, there are no con- 
straints on the model parameters from this channel. 

Finally, for the associated production with bottom 
quarks, the cross section is also suppressed due to 
a moderate tan/3 and three-body final state phase 
space. The cross section for e+ bbai at ^/s = 

200 GeV is calculated to be 1.5 x lO'^ fb for the 
model parameters in Table HI which also indicates 
no constraints from LEP. 



B. Tevatron constraints 

The lightest CP-even Higgs boson hi has approx- 
imately the same couplings to fermions as in the 
Standard Model. The ongoing searches at CDF and 
DO do not yet constrain an hi with a mass around 
115 GeV at the Tevatron [H]. 

For the lightest CP-odd Higgs boson ai, the main 
production process at the Tevatron is through asso- 
ciated production with bb and has a cross section 
proportional to tan^ P cos^ 9a (^{bb (j>sM), where cj) 
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has the same couphng to 65 as the SM Higgs bo- 
son. Ahhough the production cross section can be 
large for the model parameters that we consider, the 
existing searches at DO for the oi decays to two taus 
only constrain an ai with a mass above 90 GeV j45j . 

However, the recent searches for ai in the chan- 
nel hi — *■ oioi —^ jX^ ii~ jX^ ^~ at DO can impose 
a stringent bound on the muon-favored model pa- 
rameter space in the NMSSM. The SM background 
for two pairs of coUinear muons is very small (below 
0.02 events for 3.7 fb~^ integrated luminosity), and 
the null result imposes a constraint a{pp — > hi -\-X) ■ 
Br(/ii ^ aiai) •Br(ai ^ M+M")^ ^ 10 fb [Mj- This 
bound is roughly independent of the ai mass. The 
production cross section for a SM Higgs at the Teva- 
tron is around 0.8 pb, so we need to have Br(/ii — 
aifli) < 1.2% for Br(ai ^ ~ 100%. There- 

fore, the muon-favored model point in Table U is al- 
lowed by this constraint, although it selects a specific 
region of the NMSSM parameter space. This exper- 
imental result points to small values of A and k to 
decrease the hi branching ratio into aioi. However, 
small values of A and k also increase the branching 
ratio of 02 into bb. Since the dominant annihila- 
tion is mediated by the resonance effect with 02, 
this leads to a non-neglegible hadronic final state 
from the dark matter annihilation. 



C. The magnetic moment of the muon 

There is a new radiative contribution to the mag- 
netic moment of the muon by exchanging oi in 
the l oop diagram. Using the one-loop result from 
Ref. [431 , the new physics correction to is 



7 9 
92 mt 



32 7r2 M, 

,2 



W 
1 



tan/3)2 
—x^ dx 



x^ml/Ml) + 1 



(75) 



and is negative. The two-loop calculation will not 
change the sign of Sa^ for a mass of oi below 
1 GeV [l3|. The measured value of has a 3.4 a 
deviation (e^ data only) above the prediction of 



the standard model [48| . Requiring the new physics 
to be less than the experimental error, we arrive at 
the following constraints on the model parameters: 

cos6lAtan^ < 2.5 (BNL (5-2)^/2). (76) 

Here Ma-^ is chosen to be 800 MeV, while the con- 
straints are less stringent as one increases Ma-^ ■ 

Notice that the Sa^ may also receive significant 
contributions from other particles like smuons, that 
we chose to be heavy here to isolate the dark matter 
discussions. 



D. Constraints from Upsilon decays 

Another stringent bound on the NMSSM param- 
eter space with a light oi below 10 GeV is from 
Upsilon decays into a photon plus oi, which decays 
into a pair of taus or muons. 

For the mass range 2mr < < 9.2 GeV, the 

strongest bound is from the recent CLEO-III lim- 
its [431 on T(15') "/T~^T~. The radiative decay 
to T(IS') ^ 7ai is calculated as [50l | 



G 



B{T{1S) ^ 7ai) _ 
B{T{1S) ^ fi+ fi-) V2^< 



{cos 9 A tan/3)^ 



Ml 
Mr (IS) , 



(77) 



where T ~ 0.5 incorporates QCD and relativistic 
corrections [5l| [12] ■ The data from CLEG have 
the limit B{J{\S) 7ai) x B{ai r+r") < 
5x 10~^ at 90% C.L. for a wide range of Ma^ between 
4 GeV to 9 GeV. Using S(T(1S') ^ M+M") = 
2.48% and B{ai t+t") 0.9 from the model 
prediction, this limit is translated into a bound on 
cos Ba tan /? as 

cosQa tan/3 < 0.9 (CLEG - HI) . (78) 

A similar result is obtained in [5^. For the mass 
range 2m^ < < 1 GeV, the strongest current 

bound is coming from the light scalar searches in 
the channel T(3S') ^ 7ai at BaBar. At 90% C.L., 
BaBar imposes an upper limit S(T(3S') — s- 701) x 
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for below 



B{ai -> M+Ai") ^ 5.2 X lO"*^ | 
1 GeV. Using B{f{iS) H+ = 2.18% and 
B{ai — ^ ~ 1.0 from the model prediction, 

this limit is translated into a bound on cos 6a tan (3 
as 



cos 6*^ tan/3 < 0.4 (BaBar). 



(79) 



For the range of 2m ^ < niai < — M^^, the 
decay mode 7r+ ai is open. The branching 

ratio of this decay channel is given by [55| 



1^ ai) 



V2 

X cos^ 6a 



(tan/3 - i&ir^ PY 



r(if+ ^ all) 



= 3 X lO-'^ (tan/3 - tan^i /3)2 cos^ 6^ 



(80) 



which should be compared to the experimental val- 
ues from the HyperCP collaboration 56]: B{K^ — > 
7r+/i+^-) = 9.8 ± 1.0 ± 0.5 X IQ-*^. Therefore, the 
following constraint on the model parameter space 
is derived 

cos^A I tan ^ - tan-i ^| < 0.06 (HyperCP)(31) 

So, other than when tan /3 is very close to 1 or cos 6a 
is extremely close to zero, the mass of ai is con- 
strained to be above Mk — M.„- ~ 360 MeV. 

Finally, there are also other constraints from B- 
physics like b —^ sj oi —^ fi'^ ^~ . Since the 
oi pseudoscalar does not mediate tree-level flavor- 
changing processes, one can use the minimal flavor 
violation assumption to suppress many kinds of fla- 
vor changing processes. We have used the NMHDE- 
CAY program to check those constraints and to jus- 
tify the validity of our model points. In short, the 
most stringent constraint for Mk — M^^ < nia^ < 
1 GeV is the radiative decays of Upsilon into pho- 
tons from BaBar. The bound is cos 6 a tan/3 < 0.4. 
One can see that the model point in Table [J satisfies 
this bound. 



VIII. DISCUSSIONS AND CONCLUSIONS 

The dark matter candidate neutralino in the 
NMSSM from Table |T] is a combination of the bino. 



wino and Higgsino. Therefore, it has a good chance 
to be detected in dark matter direct detection ex- 
periments, especially from the spin-dependent elas- 
tic scattering with nucleons. Here we just report 
the values calculated using the micrOMEGAs pro- 
gram (srj . The spin- independent DM-proton and 
DM-neutron cross sections are 0.7 x 10^*^ cm^ and 
0.9 X 10"'*^ cm^, which are two orders of magni- 
tude below the current bounds from XENONIO [58j . 
For the spin-dependent one, the DM-proton and 
DM-neutron cross sections are 1.5 x 10~^^ cm^ and 
1.2 X 10~^^ cm^. The later one is only one order of 
magnitude below the current bound [595 ^^'^ 
the accessible region of the upgraded experiments. 

The NMSSM explanation of the PAMELA 
positron excess can be tested by future experiments 
at colliders and new results from cosmic ray experi- 
ments. On the collider side, it is important to mea- 
sure the masses of the neutralino and the heavier 
CP-odd particle. If their masses satisfy the relation 
~ 2m^, the large dark matter annihilation 
cross section can be confirmed. Another important 
quantity to measure is the mass of the lighter CP- 
odd particle, because the PAMELA positron excess 
prefers to have its mass below 1 GeV. It is also cru- 
cial to know the branching ratio of 02 to bb, since 
the dominant antiproton contributions are from this 
channel. 

On the cosmic ray side, we make a well-motivated 
unambiguous prediction that PAMELA will observe 
a turnover of the rising positron spectrum at around 
70 GeV. The additional contributions to the elec- 
tron+positron spectrum from dark matter annihila- 
tions are within current uncertainties but could be 
resolved in the future by Fermi LAT. The associ- 
ated gamma-ray flux from dark matter annihilation 
could also be resolved by Fermi LAT, but again due 
to astrophysical uncertainties one should be cautious 
when making a concrete prediction. 

In this paper we have explored the possibility 
of using neutralino annihilations in the NMSSM to 
explain the positron excess observed at PAMELA. 
Kinematics plays an essential role for having a vi- 
able model with a large fraction of leptons and a 
small fraction of hadrons in the annihilation final 
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state. The lighter CP-odd particle ai has a mass 
below 1 GeV and mainly decays into two muons. 
The dark matter neutralino mass is less than the 
top quark mass to forbid the otherwise dominant 
tt final state. The s-channel resonance effect with 
the heavier CP-odd particle 02 increases the dark 
matter annihilation cross section to match the nec- 
essary one for the PAMELA positron fraction spec- 
trum. This also requires that the neutralino mass is 
less than the top quark mass, to avoid smearing the 
resonance effect from a large 1 1 contribution to the 
02 width. 

We have also shown that there is discrepancy be- 
tween the NMSSM predictions and the preliminary 
gamma ray fluxes from Fermi LAT. To alleviate this 
descrepancy, one could evoke the existence of an as- 
trophysical boost factor from a nearby clump of dark 
matter. Such boost factor would also improve the 
agreement between the NMSSM prediction for the 
antiproton spectrum and the PAMELA results. In 
particular, given the fact that positrons/electrons 
are only coming from nearby sources, a clump of 
dark matter can affect positron/electron, antipro- 
ton/proton and gamma ray fluxes differently, and 
would relatively increase the positron ratio spectrum 
more than the antiproton ratio spectrum and the 
gamma ray flux spectrum. We have not utilized this 



astrophysical boost factor in our analysis. 

An NMSSM explanation of PAMELA makes three 
striking and uniquely correlated predictions: the rise 
in the PAMELA positron spectrum will turn over at 
around 70 GeV, the dark matter particle mass is 
less than the top quark mass, and a light sub-GeV 
pseudoscalar will be discovered at colliders. 
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Note: After this paper was submitted in the 
arXiv, the BaBar collaboration presented a new con- 
straint on the Upsilon radiative decays to photon 
plus two muons by combining both T(2S) and T(3S') 
data fsd]. This implies a more stringent bound on 
cos 0^ tan/? by a factor of two, see Eq. |yff[ ). This 
can be accommodated in the NMSSM by increasing 
the branching ratio of 02 to b b. This increases the 
tension for the antiproton ratio spectrum observed 
at PAMELA. However, a modest astrophysical boost 
factor on the order 0/ < 5 would loosen this tension. 
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